1.1 NVIDIA

Sumber : finance.yahoo.com

Terdapat beberapa faktor yang menjadi landasan dan membuat kami memilih saham NVIDA. Berikut adalah latar belakang yang melandasi pemilihan ini:

1. Dominasi di Industri Chip AI
   NVIDIA Corporation tidak hanya merupakan pemimpin, tetapi juga menggambarkan dominasi dalam industri chip kecerdasan buatan (AI). Sejak pendiriannya pada tahun 1993, perusahaan ini telah memainkan peran kunci dalam menghadirkan solusi grafis dan komputasi untuk berbagai industri, termasuk gaming, visualisasi profesional, pusat data, dan otomotif.
2. Inovasi Chip H200 sebagai Pemacu Utama
   Pengumuman pengembangan chip terbaru, H200, menambahkan dimensi baru pada potensi inovatif NVIDIA dalam mendukung perkembangan AI. Dengan peningkatan signifikan pada memori dan kecepatan, H200 memberikan daya ungkit yang kuat bagi aplikasi AI generatif dan model bahasa besar, menjadikannya subjek yang menarik untuk diteliti.
3. Peran Penting dalam Layanan AI Generatif
   NVIDIA tidak hanya mendominasi pasar chip AI, tetapi juga memainkan peran penting dalam layanan AI generatif, termasuk layanan ChatGPT yang populer. Sebagai penyedia teknologi di balik berbagai aplikasi AI canggih, perusahaan ini memiliki dampak yang signifikan pada perkembangan dan kemajuan kecerdasan buatan.
4. Dukungan dari Cloud Service Providers Terkemuka
   Keterlibatan NVIDIA dengan penyedia Cloud Service terkemuka seperti Amazon Web Services, Google Cloud, Microsoft Azure, dan Oracle Cloud Infrastructure menunjukkan pengakuan industri terhadap kualitas dan relevansi produk-produknya. Ini menambah nilai pada data saham NVIDIA sebagai subjek penelitian.
5. Performa Saham yang Mengesankan
   Kinerja saham NVIDIA yang mengesankan, terutama dengan kenaikan lebih dari tiga kali lipat dalam nilai selama tahun ini, menarik perhatian sebagai indikator potensial untuk menganalisis dampak inovasi teknologinya terhadap nilai pasar saham.

Dengan kombinasi faktor-faktor ini, pemilihan data saham NVIDIA bukan hanya memungkinkan analisis yang mendalam tentang hubungan antara inovasi teknologi dan pergerakan harga saham, tetapi juga memberikan wawasan yang berharga tentang bagaimana NVIDIA terus memainkan peran krusial dalam perkembangan industri kecerdasan buatan.

1.2 Dataset

Dataset yang saya gunakan merupakan koleksi data harga saham historis periode 1 Januari 2018 hingga 10 November 2023 dari 6 saham yakni saham saingan NVDA seperti AMD, ARM, AVGO, TSM, dan INTC.

Dataset ini memilki data :

1. Open: yakni Harga saham pada awal periode perdagangan tertentu. Ini adalah harga saham pertama pada hari perdagangan tersebut.
2. High: Harga tertinggi yang saham capai selama periode perdagangan tersebut. Ini mencerminkan harga tertinggi yang pembeli bersedia bayar selama hari tersebut.
3. Low: Harga terendah yang saham capai selama periode perdagangan tersebut. Ini mencerminkan harga terendah yang penjual bersedia terima selama hari tersebut.
4. Close: Harga saham pada akhir periode perdagangan tertentu. Ini adalah harga saham terakhir pada hari perdagangan tersebut.
5. Adj Close (Adjusted Close): Harga penutup yang telah disesuaikan untuk memperhitungkan perubahan seperti pembagian saham atau dividen. Ini adalah harga penutup yang paling relevan untuk analisis jangka panjang, karena mencerminkan harga saham yang sebenarnya setelah penyesuaian.
6. Volume: Volume perdagangan saham selama periode tertentu. Ini mencerminkan jumlah saham yang diperdagangkan selama hari perdagangan tersebut.

Kami akan menggunakan peubah Adj Close (Adjusted Close), Karena sesuai dengan penjelasan diatas, peubah Adj Close adalah yang paling sesuai untuk dianalisis dibandingkan peubah lainnya.

1.3 Tujuan

Tujuan dari praktikum ini adalah untuk menganalisis pola perkiraan pergerakan harga tujuh saham teknologi terkemuka dengan harapan dapat memberikan rekomendasi kepada pembaca mengenai saham mana yang sebaiknya dipertimbangkan untuk dibeli atau diinvestasikan secara signifikan di antara tujuh perusahaan teknologi besar tersebut.

2.1 Plot Time Series

install_load('ggplot2','extrafont')
# font_import(); loadfonts() #Run ini sekali aja
theme.ts <- list(
  theme(legend.position = "none",
        axis.text.x = element_text(hjust = 1, 
                                   margin = margin(b = 10, t=20)),
        axis.text.y = element_text(vjust = 0.5, face = "bold", 
                                   margin = margin(l = 20, r = 20)),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        text = element_text(size = 30),
        plot.subtitle = element_text(hjust = 0.5),
        panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill='transparent', color=NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(linewidth = 1, colour = "black"))
        )
theme.ts1 <- list(
  theme(legend.position = "none",
        axis.text.x = element_text(hjust = 1, 
                                   margin = margin(b = 10, t=20)),
        axis.text.y = element_text(vjust = 0.5, face = "bold", 
                                   margin = margin(l = 50, r = 20)),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        text = element_text(size = 30),
        plot.subtitle = element_text(hjust = 0.5),
        panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill='transparent', color=NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(linewidth = 1, colour = "black"))
        )

install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
  geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
  #Color
  scale_color_manual(values = c(AMD=col.amd[1],   ARM=col.arm[1],
                                AVGO=col.avgo[1], NVDA=col.nvda[1],
                                TSM=col.tsm[1],   INTC=col.intc[1] )) +
  scale_alpha_manual(values = c("AMD" = .25, "ARM" = .25, "AVGO" = .25, 
                                "NVDA" = 1, "TSM" = .25, "INTC" = .25)) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series NVDA",
       subtitle = "Seperti apa sih pola deret waktu saham NVDA dan saingannya?\n") +
  # Label / legend
  geom_text_repel(
    data=data[data$Date == max(data$Date),], #Posisi di ujung data
    aes(color = Name, label = Name), #Warna garis & label saham
    size = 8, #Ukuran text
    nudge_x = 80, #Posisi Text (kanan 50)
    hjust = 0, #Ujung
    segment.size = 1,               #Ukuran garis
    segment.alpha = .75,             #transparasi garis
    segment.linetype = "dotted",    #Time garis
    box.padding = .4, #Biar label saham nggak dempetan
    segment.curvature = -0.1, #biar garis mulus
    segment.ncp = 8, 
    segment.angle = 60 
  ) +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(data$Date)), 
               as.Date(max(data$Date)) + 120)
    #Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
    annotate( #Buat nandain batas data
    "text", x = as.Date(max(data$Date)), y = -35, 
    label = max(data$Date), size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data$Date) )), 
             linetype = "dotted", color = "red")
chart

#Export Chart
ggsave("01_Time Series AAANTI.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

data2 <- data %>% 
  filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10")) 

install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data2, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
  geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
  #Color
  scale_color_manual(values = c(AMD=col.amd[1],   ARM=col.arm[1],
                                AVGO=col.avgo[1], NVDA=col.nvda[1],
                                TSM=col.tsm[1],   INTC=col.intc[1] )) +
  scale_alpha_manual(values = c("AMD" = .35, "ARM" = .35, "AVGO" = .35, 
                                "NVDA" = 1, "TSM" = .35, "INTC" = .35)) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series NVDA 2021-2023",
       subtitle = "Seperti apa sih pola deret waktu saham NVDA dan saingannya?\n") +
  # Label / legend
  geom_text_repel(
    data=data2[data2$Date == max(data2$Date),], #Posisi di ujung data
    aes(color = Name, label = Name), #Warna garis & label saham
    size = 8, #Ukuran text
    nudge_x = 20, #Posisi Text (kanan 50)
    hjust = 0, #Ujung
    segment.size = 1,               #Ukuran garis
    segment.alpha = .75,             #transparasi garis
    segment.linetype = "dotted",    #Time garis
    box.padding = .4, #Biar label saham nggak dempetan
    segment.curvature = -0.1, #biar garis mulus
    segment.ncp = 8, 
    segment.angle = 60 
  ) +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(data2$Date)), 
               as.Date(max(data2$Date)) + 60)
    #Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
    annotate( #Buat nandain batas data
    "text", x = as.Date(max(data2$Date)), y = -35, 
    label = max(data2$Date), size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data2$Date) )), 
             linetype = "dotted", color = "red")
chart

#Export Chart
ggsave("01_Time Series AAANTI_2022-2023.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

nvda <- data2 %>%
  filter(Name == "NVDA")  # Filter data saham Amazon tahun 2022 ke atas

rownames(nvda) <- NULL
str(nvda)

## 'data.frame':    1044 obs. of  3 variables:
##  $ Name     : chr  "NVDA" "NVDA" "NVDA" "NVDA" ...
##  $ Date     : Date, format: "2021-01-01" "2021-01-02" ...
##  $ Adj Close: num  130 131 131 131 134 ...

datatable(nvda, filter = 'top', 
          options = list(pageLength = 5))

nvda.ts <- ts(nvda[,3])

summary(nvda.ts)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   112.2   158.1   206.4   237.1   278.5   493.5

install_load('ggtext')
min_value <- min(nvda$`Adj Close`)
min_date <- nvda$Date[which.min(nvda$`Adj Close`)]
min_perc <- (which.min(nvda$`Adj Close`) / nrow(nvda)) * 100

max_value <- max(nvda$`Adj Close`)
max_date <- nvda$Date[which.max(nvda$`Adj Close`)]
max_perc <- (which.max(nvda$`Adj Close`) / nrow(nvda)) * 100

chart <-
ggplot(nvda, aes(x=Date, y=`Adj Close`)) + 
  geom_line(aes(color=Name), linewidth=2) +
  scale_color_manual(values = col.nvda[1]) +
  labs(x = "\nPeriode (Tahun)", y='Saham Harga penutup',
       title = "Time Series Saham NVIDIA",
       subtitle = "Seperti apa sih pola deret waktu saham NVIDIA?\n") +
  theme(legend.position = "none") +
  theme.ts1 + 
  
#Ekstras

  #Titik terendah
  geom_segment(aes(x = min_date, 
                   xend = min_date, 
                   y = min(nvda$`Adj Close`)* 1.01, 
                   yend = max(nvda$`Adj Close`)*40/100), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = min_date, 
                      y = max(nvda$`Adj Close`)*40/100, 
                      label = paste0("Titik Terendah") ),
    aes(x, y, label = label), size = 7, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  ) +
  geom_text(aes(x = min_date-1*35, y = max(`Adj Close`)*40/100, label = 
                  paste0(min_date)), 
            vjust = -1.5, hjust = 0, size = 7, color = "grey30") +

  #Titik tertinggi
  geom_segment(aes(x = max_date, 
                   xend = max_date, 
                   y = max(nvda$`Adj Close`)* .99, 
                   yend = max(nvda$`Adj Close`)*70/100), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = max_date, 
                      y = max(nvda$`Adj Close`)*70/100, 
                      label = paste0("Titik Tertinggi") ),
    aes(x, y, label = label), size = 7, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) + 
  geom_text(aes(x = max_date-1*35, y = max(`Adj Close`)*60/100, label = 
                  paste0(max_date)), 
            vjust = -1.5, hjust = 0, size = 7, color = "grey30") +
  annotate( #Buat nandain batas data
    "text", x = as.Date(max(data2$Date)) -32, y = 50, 
    label = max(data2$Date), size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data2$Date) )), 
             linetype = "dotted", color = "red", linewidth=1.5) 

chart

#Export Chart
ggsave("02_Time Series NVIDIA.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

install_load("forecast")

## Warning: package 'forecast' was built under R version 4.2.3

ggseasonplot(nvda.ts)

## Error in ggseasonplot(nvda.ts): Data are not seasonal

2.2 Data Train vs Test

#membagi 80% data latih (training) dan 20% data uji (testing)
train_nvda <- nvda[1: round(nrow(nvda) *78/100),]
test_nvda <- nvda[round(nrow(nvda) *78/100  +1): nrow(nvda),]
train_nvda.ts <- ts(train_nvda[,3])
test_nvda.ts <- ts(test_nvda[,3])

chart <-
ggplot() + 
  #Label Data Asli 
  annotate( "rect", alpha=.1, fill="#4EC2C1",
            xmin=as.Date(min(data2$Date)), 
            xmax=as.Date(data2$Date[nrow(train_nvda)]),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Latih", size=10) + 
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(train_nvda)," Hari (", 
                   round(nrow(train_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) + 
  
  #Label Data Ramal
  annotate( "rect", alpha=.1, fill="violetred",
            xmin=as.Date(data2$Date[nrow(train_nvda)]), 
            xmax=as.Date(max(data2$Date)),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Uji", size=10) +
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(test_nvda)," Hari (", 
                   round(nrow(test_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) +
  
  geom_vline( #Buat garis batas data
    xintercept = as.Date(data2$Date[nrow(train_nvda)]) , 
             linetype = "dotted", color = "black", linewidth = 2) +

#NVDA
  geom_line(data = train_nvda, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "NVDA")) +
  geom_line(data = test_nvda, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "test_nvda")) +
  geom_point(data = tail(train_nvda, 1), alpha = .5, 
             aes(x = Date, y = `Adj Close`), stroke=2,
             size = 15, shape = 21, color = "black", fill="violetred") +
  scale_colour_manual(values = c("NVDA"=col.nvda[1], "test_nvda"=col.nvda[2])) +
    theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series Saham NVDA",
       subtitle = "Pembagian Data Latih dan data Uji\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +

#Bagi Data Train & Test
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)), 
                   y = max(nvda$`Adj Close`)* .58, 
                   yend = max(nvda$`Adj Close`)*.7), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)), 
                      y = max(nvda$`Adj Close`)*.7, 
                      label = max(train_nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +

  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)), 
                   y = max(nvda$`Adj Close`)* .95, 
                   yend = max(nvda$`Adj Close`)*1.05), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)), 
                      y = max(nvda$`Adj Close`)*1.05, 
                      label = max(data2$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  )
  
chart

#Export Chart
ggsave("02_TS_NVDA_train-test.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Berdasarkan plot data deret waktu pada data latih (\(78\%\) dari data asli), terlihat bahwa data menunjukkan tren naik dan tren turun. Ini mengisyaratkan bahwa data latih tidak memenuhi kriteria stasioneritas dalam rataan maupun ragam. Di sisi lain, dalam plot data uji (\(22\%\) dari data asli), terlihat adanya tren yang melonjak naik dan kurangnya nilai tengah yang stabil. Ini juga menunjukkan bahwa data uji tidak stasioner dalam rataan.

3.1 Uji Stasioneritas

install_load('tsibble','tseries')
acf(train_nvda.ts, main="", col=col.nvda[1], lwd=2)
mtext("NVDA", side=3, line=1, cex=2, font=2)

tseries::adf.test(train_nvda.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_nvda.ts
## Dickey-Fuller = -1.6692, Lag order = 9, p-value = 0.7183
## alternative hypothesis: stationary

if (tseries::adf.test(train_nvda.ts)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
  cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}

## Karena p-value > 0.05, Maka Rataan Tidak Stasioner

install_load('MASS')
index <- seq(1:nrow(train_nvda)) #Beda titik potong
bc =boxcox(train_nvda.ts~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main = 2)

lambda_nvda <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK

lsk <- data.frame(
  Lambda = c(lambda_nvda ),
  Batas.Bawah = c( min(sk_nvda)  ),
  Batas.Atas = c(max(sk_nvda) )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk$Keterangan <- apply(lsk, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk) <- c("NVDA")

datatable(lsk, filter = "none")

3.2 Penanganan Ketidakstasioneran Data

lamb <- c(-2, -1, -.5, 0, .5, 1, 2)
# Fungsi untuk mengonversi lambda ke transformasi
trans <- function(lambda) {
  lambda <- lamb[sapply(lambda, function(y) which.min(abs(y - lamb))) ]
  case_when(
    lambda == -2   ~ "1/x^2",
    lambda == -1   ~ "1/x",
    lambda == -0.5 ~ "1/sqrt(x)",
    lambda == 0    ~ "log(x)",
    lambda == 1    ~ "x",
    lambda == 0.5  ~ "sqrt(x)",
    lambda == 2    ~ "x^2",
    TRUE           ~ NA_character_
  )
}

# Fungsi untuk mengonversi lambda ke transformasi balik
trans_balik <- function(lambda) {
  lambda <- lamb[sapply(lambda, function(y) which.min(abs(y - lamb))) ]
  case_when(
    lambda == -2   ~ "1/sqrt(x)",
    lambda == -1   ~ "1/x",
    lambda == -0.5 ~ "1/x^2",
    lambda == 0    ~ "exp(x)",
    lambda == 1    ~ "x",
    lambda == 0.5  ~ "x^2",
    lambda == 2    ~ "sqrt(x)",
    TRUE           ~ NA_character_
  )
}

# Menambahkan kolom Transformasi dan Trans Balik
tran <- lsk %>%
  mutate(
    Transformasi = trans(Lambda),
    `Transformasi Balik` = trans_balik(Lambda)
  ) %>% dplyr::select(Lambda, Transformasi, `Transformasi Balik`)

datatable(tran, filter = 'none')

#Transformasi
train_nvda.ts.new <- 1/sqrt(train_nvda.ts)

#Plot BoxCox
index <- seq(1:nrow(train_nvda))
bc <-boxcox(train_nvda.ts.new ~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main=2)

lambda_nvda.new <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.new <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK

lsk.new <- data.frame(
  Lambda = c(lambda_nvda.new ),
  Batas.Bawah = c( min(sk_nvda.new)  ),
  Batas.Atas = c( max(sk_nvda.new)  )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk.new$Keterangan <- apply(lsk.new, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk.new) <- c("NVDA")

datatable(lsk.new, filter = 'none')

ragam <- data.frame(
  Date = train_nvda %>% dplyr::select(Date)  ,
  `Asli` = c(as.vector(train_nvda.ts) ),
  `Trans` = c(as.vector(train_nvda.ts.new) ) *3000,
  Name = rep(c("NVDA"), 
              each = nrow(train_nvda) )
)
colnames(ragam) <- c("Date", "Asli", "Trans", "Name")

chart <- 
ggplot() +
# Akhir data
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(train_nvda$Date) )), 
             linetype = "dotted", color = "red") +

#Time Series
  #Data Asli
  geom_line(data = ragam, linewidth=2,
            aes(x = Date, y = `Asli`, col = "Asli")) +
  #Data Transformasi
  geom_line(data = ragam, linewidth=2,
            aes(x = Date, y = `Trans`, col = "Trans")) +

  scale_colour_manual(values = c("Asli"=col.nvda[2], "Trans"=col.nvda[1])) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Plot Transformasi NVDA",
       subtitle = "Perbandingan plot Transformasi dengan Plot Asli\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(ragam$Date)), 
               as.Date(max(ragam$Date)) + 80)
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
  ylim(min(ragam$Asli), max(ragam$Asli)) +

#Data asli
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date))+60, 
                   y = train_nvda$`Adj Close`[nrow(train_nvda)], 
                   yend = train_nvda$`Adj Close`[nrow(train_nvda)]), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date))+60, 
                      y = train_nvda$`Adj Close`[nrow(train_nvda)], 
                      label = "Data Asli" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +

# Data Transformasi
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)) + 60, 
                   y = ragam$Trans[nrow(ragam)], 
                   yend = ragam$Trans[nrow(ragam)] ), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)) +60, 
                      y = ragam$Trans[nrow(ragam)], 
                      label = "Transformasi" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  ) +
  
#Tanggal akhir
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)), 
                   y = -Inf, 
                   yend = max(train_nvda$`Adj Close`)*.15), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)), 
                      y = max(train_nvda$`Adj Close`)*.15, 
                      label = max(train_nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  )

chart

#Export Chart
ggsave("03_Stas_4_Transform.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

#Diff
train_nvda.diff <- diff(train_nvda.ts.new, differences = 1) 

#Plot
plot.ts(train_nvda.diff, lty=1, xlab="Periode (Hari)", 
        col = col.nvda[1], lwd = 2,
        main="Plot Differencing", cex.main = 2)

3.3 Uji Ulang

Akan dilakukan uji ulang untuk mengecek kestasioneran data dalam rataan.

stats::acf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)

tseries::adf.test(train_nvda.diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_nvda.diff
## Dickey-Fuller = -8.3563, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

if (tseries::adf.test(train_nvda.diff)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
  cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}

## Karena p-value < 0.05, Maka Rataan Stasioner

4.1 Model Tentatif

Identifikasi model tentatif dapat dilakukan dengan melihat pada plot ACF, PACF dan Tabel EACF. Mengutip dari postingan pada laman medium.com Panduan dasar untuk menginterpretasikan plot ACF dan PACF adalah sebagai berikut :

1. Cari pola tail off pada ACF atau PACF.
2. Jika tail off terjadi pada ACF → model AR → cut off pada PACF akan memberikan nilai p untuk AR(p).
3. Jika tail off terjadi pada PACF → model MA → cut off pada ACF akan memberikan nilai q untuk MA(q).
4. Jika terjadi tail off pada kedua ACF dan PACF → model ARMA.

par(mfrow=c(1,2))
stats::acf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)

stats::pacf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)

acf.pacf <- function(dt){ #Fungsi Deteksi Cut off plot ACF & PACF
  N <- acf(dt, plot = F)[["n.used"]] + 1 #Banyaknya Data
  #Garis Signifikan
  sl <- c(-1,1)*(exp(2*1.96/sqrt(N-3))-1)/(exp(2*1.96/sqrt(N-3))+1)
  acf <- data.frame( lag = acf(dt, plot=F)[["lag"]], 
                     acf = acf(dt, plot=F)[["acf"]] ) 
  
    cut.off1 <- acf %>%
        filter(acf <= min(sl) | acf >= max(sl)) %>%
        dplyr::select(lag); colnames(cut.off1) <- "acf"
  
  pacf <- data.frame(lag = pacf(dt, plot=F)[["lag"]], 
                     acf = pacf(dt, plot=F)[["acf"]] ) 
  
    cut.off2 <- pacf %>%
        filter(acf <= min(sl) | acf >= max(sl)) %>%
        dplyr::select(lag); colnames(cut.off2) <- "pacf"
    
  cut.off <- merge(cut.off1, cut.off2, by = 0, all = TRUE)[-1]
  return(cut.off)
}

data_frames <- list(
  train_nvda.diff
)

name <- "NVDA"

acf.pacf.tab <- data.frame()
for(i in 1:length(name)){
  dex <- data_frames[[i]] %>% acf.pacf() %>% 
  rename(!!paste0(name[i], ".a") := acf, 
         !!paste0(name[i], ".p") := pacf)
  
    acf.pacf.tab <- merge(acf.pacf.tab, dex, by=0, all=TRUE)[-1]
}

datatable(acf.pacf.tab, filter = 'none')

install_load('TSA', 'forecast')

# Fungsi untuk menghitung model ARIMA dan menganalisis parameter
model_tentatif <- function(data, p_max, d, q_max, alpha=0.05, drift=FALSE) {
  best_model <- NULL
  best_aic <- Inf
  eacf_result <- eacf(data) #Untuk Ektraksi & identifikasi model
  models <- data.frame(Model = character(0), 
                       AIC = numeric(0), 
                       Signif = character(0), 
                       Keterangan = character(0))
  
  for (p in 0:p_max) {
    overfit <- 0
    for (q in 1:q_max) {
      #Pola Matriks segitiga bawah
      if (!is.na(eacf_result$symbol[p + 1, q ]) && 
          !is.na(eacf_result$symbol[p + 1, q + 1]) && 
          !is.na(eacf_result$symbol[p + 2, q + 1])) {
        if (eacf_result$symbol[p + 1, q ] == "o" && 
            eacf_result$symbol[p + 1, q + 1] == "o" && 
            eacf_result$symbol[p + 2, q + 1] == "o") {
      
          model <- Arima(data, order=c(p,d,q), method="ML", 
                         include.drift=drift)
          aic <- AIC(model)
          
          # Mendapatkan nilai coef dari model
          coeftest_result <- lmtest::coeftest(model)
          
          # jika lebih kecil dari alpha, maka signifikan
          significant_params <- 
            rownames(coeftest_result)[coeftest_result[, "Pr(>|z|)"] < alpha]  
          
          # jika lebih besar dari alpha, maka tidak signifikan
          non_significant_params <- 
            rownames(coeftest_result)[coeftest_result[, "Pr(>|z|)"] > alpha]  
          
          # Keterangan signifikansi
          if (length(significant_params) == 0) {
            keterangan <- "Semua parameter tidak signifikan"
          } else if (length(significant_params) == nrow(coeftest_result)) {
            keterangan <- "Semua parameter signifikan"
          } else {
            keterangan <- paste("Parameter yang tidak signifikan adalah", 
                                paste(non_significant_params, collapse = ", "))
          }
          
          models <- rbind(models, 
                    data.frame(Model = paste("ARIMA(", p, ",", d, ",", q, ")", 
                                             sep = ""), 
                               AIC = aic, 
                               Signif = paste(significant_params, 
                                              collapse = ", "), 
                               Keterangan = keterangan))
          
          #Identifikasi Best Model
            if (keterangan == "Semua parameter signifikan" && 
                !any(is.na(significant_params))) {
              if (aic < best_aic) {
                best_model <- model
                best_aic <- aic
            }
          }
          #jika model ditemukan pada baris p. Maka lanjut ke p+1
          if(overfit >= 1) break
          overfit = overfit + 1
        }
      }
    }
  }
  
  cat("\nModel ARIMA dengan AIC terkecil:\n")
  print(best_model)

  return(models)
}

model.tentaif_nvda <- 
  model_tentatif(train_nvda.diff, p_max = 6, d = 1, q_max = 12)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o x o o o o o o  o  o  o 
## 6 x x x x x x o o o o o  o  o  o 
## 7 x x o x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## Series: data 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.9933
## s.e.   0.0078
## 
## sigma^2 = 9.868e-07:  log likelihood = 4460.68
## AIC=-8917.36   AICc=-8917.34   BIC=-8907.96

datatable(model.tentaif_nvda, filter = 'top', 
          options = list(pageLength = 7))

model_nvda.da <- Arima(train_nvda.diff, order=c(0,1,1), method="ML")

4.2 Model ARIMA dengan Drift

Tren menentukan bagaimana perubahan deret waktu dari waktu ke waktu. Dalam model ARIMA, ini dapat dimodelkan menggunakan persamaan dasar berikut “a” disebut “intercept” dan “b” disebut “drift”. “Drift” adalah hanyalah kemiringan garis lurus. Selanjutnya akan dilakukan kombinasi model ARIMA dengan Intercept dan Drift.

drift1_nvda <- 
  model_tentatif(train_nvda.ts, p_max = 6, d = 1, q_max = 12)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 o o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x x o o o o o o o o  o  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x x o x x o o o o o o  o  o  o 
## 7 x x x x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## Series: data 
## ARIMA(4,1,4) 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1     ma2      ma3     ma4
##       0.4960  -0.7914  0.5091  -0.9828  -0.4804  0.7803  -0.5046  0.9699
## s.e.  0.0144   0.0120  0.0094   0.0172   0.0236  0.0193   0.0172  0.0240
## 
## sigma^2 = 29.16:  log likelihood = -2522.64
## AIC=5063.29   AICc=5063.51   BIC=5105.59

datatable(drift1_nvda, filter = 'top', 
          options = list(pageLength = 7))

model_nvda.drift1 <- Arima(train_nvda.ts, order=c(4,1,4), method="ML")

drift2_nvda <- 
  model_tentatif(train_nvda.diff, p_max = 6, d = 0, q_max = 12)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o x o o o o o o  o  o  o 
## 6 x x x x x x o o o o o  o  o  o 
## 7 x x o x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## NULL

datatable(drift2_nvda, filter = 'top', 
          options = list(pageLength = 7))

model_nvda.drift2 <- NA

drift3_nvda <- 
  model_tentatif(train_nvda.ts, p_max = 6, d = 1, q_max = 12, drift=TRUE)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 o o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x x o o o o o o o o  o  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x x o x x o o o o o o  o  o  o 
## 7 x x x x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## NULL

datatable(drift3_nvda, filter = 'top', 
          options = list(pageLength = 7))

model_nvda.drift3 <- NA

drift4_nvda <- 
  model_tentatif(train_nvda.diff, p_max = 6, d = 1, q_max = 12, drift=TRUE)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o x o o o o o o  o  o  o 
## 6 x x x x x x o o o o o  o  o  o 
## 7 x x o x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## NULL

datatable(drift4_nvda, filter = 'top', 
          options = list(pageLength = 7))

model_nvda.drift4 <- NA

4.3 Model Terbaik

model_terbaik <- function(model, rowname, mods) {
  if (!all(is.na(model))) {
    model_df <- data.frame(
      Model = paste("ARIMA(", 
                     paste(as.character(model[["call"]][["order"]][-1]), 
                           collapse = ","), ")", sep = ""),
      AIC = model[["aicc"]],
      row.names = rowname
    )
    model_df$Keterangan <- 
      mods[which(mods$Model == 
                     model_df["Model"][rowname,]), "Keterangan"]
  return(model_df)
    
  } else {
    return(data.frame(
      Model = NA_character_,
      AIC = NA_real_,
      Keterangan = NA_character_,
      row.names = rowname
    ))
  }
}

best_model <- rbind(
  model_terbaik(model_nvda.da, 'NVDA', model.tentaif_nvda),
  model_terbaik(model_nvda.drift1, 'NVDA', drift1_nvda),
  model_terbaik(model_nvda.drift2, 'NVDA', drift2_nvda),
  model_terbaik(model_nvda.drift3, 'NVDA', drift3_nvda),
  model_terbaik(model_nvda.drift4, 'NVDA', drift4_nvda)
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

datatable(best_model, filter = 'top', 
          options = list(pageLength = 7))

Terlihat bahwa model paling terbaik adalah model tentatif tanpa kombinasi drift. Sangat disayangkan, dari 4 kombinasi drift yang ada, hanya kombinasi 1 yang memiliki kategori model terbaik. Sehingga model yang akan digunakan adalah model_nvda.da.

4.4 Overfitting

Overfitting dilakukan dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model terpilih ARIMA(0,1,1) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(1,1,1) dan ARIMA(0,1,2).

Perlu diingat bahwa sebelumnya saya menjelaskan bahwa fungsi model_tentatif yang kami buat itu sudah sekaligus overfitting, dan itu juga sudah membandingkan mana model terbaiknya. Sehingga sebenarnya tidak perlu dilakukan overfitting. Namun saya akan tunjukkan kembali.

datatable(model.tentaif_nvda, filter = 'top', 
          options = list(pageLength = 2))

Terlihat bahwa hanya terdapat model ARIMA(0,1,2) saja. Itu dikarenakan model ARIMA(1,1,1) tidak tersedia di Tabel EACF, sehingga perlu dilakukan manual.

Arima(train_nvda.diff, order=c(1,1,1), method="ML") %>% summary() 

## Series: train_nvda.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.0040  -0.9934
## s.e.  0.0357   0.0081
## 
## sigma^2 = 9.879e-07:  log likelihood = 4460.68
## AIC=-8915.37   AICc=-8915.34   BIC=-8901.27
## 
## Training set error measures:
##                         ME         RMSE          MAE MPE MAPE      MASE
## Training set -6.096225e-06 0.0009921142 0.0006814358 Inf  Inf 0.7812342
##                      ACF1
## Training set -0.001297599

memiliki AIC yang lebih besar dan tidak semua parameternya signifikan. Sehingga tetap model terbaiknya adalah ARIMA(0,1,1).

arima.best_model_nvda <- model_nvda.da

arma.order <- c(
  arima.best_model_nvda[["call"]][["order"]][[2]],
  arima.best_model_nvda[["call"]][["order"]][[4]]
)

5.1 Eksplorasi Sisaan

#Eksplorasi 
sisaan_nvda.da <- model_nvda.da$residuals 
par(mfrow=c(2,2)) 
# QQ-Plot
qqnorm(sisaan_nvda.da) 
qqline(sisaan_nvda.da, col = "dodgerblue3", lwd = 2.5) 
# Plot sisaan
plot(c(1:length(sisaan_nvda.da)), sisaan_nvda.da, 
     main = "Plot Sisaan Model ARIMA",
     xlab = "Periode",
     ylab = "Nilai Sisaan")
abline(h = 0, col = "firebrick2", lty = 2, lwd=2.5)
#Plot ACF dan PACF
stats::acf(sisaan_nvda.da, ylim=c(-.1,.1)) 
stats::pacf(sisaan_nvda.da, ylim=c(-.1,.1)) 

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan menyebar tidak normal ditandai dengan tititk-titiknya cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Plot ACF dan PACF sisaan ARIMA(0,1,1) juga tidak signifikan pada 20 lag awal yang menandakan saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

6.1 ARCH

ARCH merupakan singkatan dari Auto-Regressive Conditional Teteroskedasticity.

1. AutoRegressive (AR): Ini merujuk pada ide bahwa variansi (volatilitas) suatu variabel bergantung pada nilai-nilai sebelumnya dalam deret waktu. Dengan kata lain, variabilitas saat ini dipengaruhi oleh nilai-nilai sebelumnya.
2. Conditional (Conditional): Ini menunjukkan bahwa perhitungan volatilitas dilakukan secara kondisional terhadap informasi sebelumnya dalam deret waktu. Dengan kata lain, variabilitas diukur dengan mempertimbangkan informasi kondisional terkini.
3. Heteroskedasticity: Ini merujuk pada sifat ketidakseragaman (Heterogen/Tidak Homogen) dalam variabilitas. Dalam konteks ARCH, heteroskedasticity menunjukkan bahwa variabilitas suatu variabel tidak konstan, tetapi bervariasi sepanjang waktu berdasarkan informasi kondisional.

Jadi, AutoRegressive Conditional Heteroskedasticity (ARCH) adalah model statistik yang digunakan untuk mengukur dan memodelkan volatilitas dalam deret waktu finansial atau ekonometrik. Model ini mengasumsikan bahwa variabilitas harga atau hasil suatu variabel dapat bervariasi seiring waktu dan dipengaruhi oleh nilai-nilai sebelumnya dalam deret waktu.

Volatilitas sendiri merujuk pada tingkat fluktuasi atau variasi harga atau hasil suatu variabel dalam periode waktu tertentu. Volatilitas tinggi menunjukkan fluktuasi yang besar dalam nilai-nilai variabel tersebut, sementara volatilitas rendah menunjukkan fluktuasi yang lebih terbatas. Dalam konteks finansial, volatilitas sering kali dianggap sebagai ukuran risiko, karena fluktuasi harga yang besar dapat menunjukkan ketidakpastian atau risiko yang tinggi.

Model ARCH adalah model yang populer karena spesifikasi variansnya dapat menangkap fitur umum dari deret waktu variabel keuangan. Model ini berguna untuk memodelkan volatilitas dan khususnya perubahan volatilitas seiring waktu. Mari kita lihat kembali pada plot NVDA untuk mengevaluasi volatilitasnya.

#Plot
plot.ts(train_nvda.diff, lty=1, xlab="Periode (Hari)", 
        col = col.nvda[1], lwd = 2,
        main="NVDA", cex.main = 2)

Seperti yang dapat Anda lihat, nilai dari deret ini berubah dengan cepat dari periode ke periode dalam suatu cara yang tampaknya tidak dapat diprediksi, yang dapat kita sebut sebagai volatil. Selain itu, ada periode-periode ketika perubahan besar diikuti oleh perubahan besar berikutnya dan periode di mana perubahan kecil diikuti oleh perubahan kecil berikutnya, merupakan bukti lain dari volatilitas.

install_load('dynlm')
# Langkah 1: Mengestimasi persamaan rata-rata r = beta + error
byd.mean <- dynlm(train_nvda.diff ~ 1)

# Langkah 2: Mengambil residual dari model sebelumnya dan mengkuadratkannya
ehatsq <- ts(resid(byd.mean)^2)

# Langkah 3: Regresikan residual kuadrat pada residual kuadrat yang tertinggal satu lag
byd.arch <- dynlm(ehatsq ~ L(ehatsq), data = ehatsq)

summary(byd.arch)

## 
## Time series regression with "ts" data:
## Start = 2, End = 813
## 
## Call:
## dynlm(formula = ehatsq ~ L(ehatsq), data = ehatsq)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.235e-06 -9.005e-07 -7.837e-07 -2.000e-07  2.875e-05 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 9.114e-07  8.638e-08  10.552   <2e-16 ***
## L(ehatsq)   7.268e-02  3.504e-02   2.074   0.0384 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.257e-06 on 810 degrees of freedom
## Multiple R-squared:  0.005282,   Adjusted R-squared:  0.004054 
## F-statistic: 4.301 on 1 and 810 DF,  p-value: 0.0384

install_load('FinTS')
lag <- 12
arch.effect <- data.frame(p.value = numeric(0), 
                          Keterangan = character(0))
for(i in 1:lag){
  p.value <- ArchTest(sisaan_nvda.da, lags=i, demean=TRUE)[["p.value"]]
  arch.effect <- rbind(arch.effect, 
      data.frame(p.value = ifelse(round(p.value, 2) == 0, 
                                  sprintf("%.2e", p.value),
                                  round(p.value, 2) ),
                 Keterangan = ifelse(p.value < 0.05, 
                                     "Signifikan", "Tidak signifikan")
                 ) )
  rownames(arch.effect)[i] <- paste0("Lag ", i)
}
datatable(arch.effect, filter = 'top', 
          options = list(pageLength = 6))

install_load('rugarch')

create_spec <- function(i, j, arma_order, mean=mean) {
  return(ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(i, j)),
    mean.model = list(armaOrder = arma_order, 
              #model tidak harus memasukkan komponen rata-rata (intercept)
                      include.mean = mean)
  ))
}


arch.garch <- function(data, pp, qq = 0, arma.order = c(0, 0), mean=FALSE) {
  arch.model <- data.frame(Model = character(0),
                            Signif = numeric(0),
                            AIC = numeric(0),
                            Sig.par = character(0),
                            Non_Sig.par = character(0))
  
  best_model <- NULL
  best_sig <- -Inf
  best_aic <- Inf
  
  if(qq == 0){
    n <- 0
  } else{
    n <- 1
  }
  for (i in 1:pp) {
    for (j in n:qq) {
        # GARCH model
        formula <- as.formula(paste("~ garch(", i, ",", j, ")", sep = ""))
        archSpec <- create_spec(i, j, arma.order, mean=mean)
  
      archFit <- ugarchfit(spec = archSpec, data = data, solver = "hybrid")
  
      if (anyNA(coef(archFit)) || anyNA(archFit@fit[["tval"]])) {
        break
      }
  
      aic <- infocriteria(archFit)[[1]] # AIC
      p_values <- archFit@fit$robust.matcoef[, 4] # P-value
      
      # persentase banyaknya param yang signif
      sig <- sum(p_values < 0.05) / length(p_values)
  
      sig_param <- names(p_values[p_values < 0.05])
      non_sig_param <- names(p_values[p_values >= 0.05])
  
      model <- if (qq == 0) {
        paste0("ARCH(", i, ")")
      } else {
        paste0("GARCH(", i, ",", j, ")")
      }
  
      arch.model <- rbind(
        arch.model,
        data.frame(
          Model = model,
          Signif = sig,
          AIC = aic,
          Sig.par = paste(sig_param, collapse = ", "),
          Non_Sig.par = paste(non_sig_param, collapse = ", ")
        )
      )
  
      # Identifikasi best model
      if (sig > best_sig) {
        if (aic < best_aic) {
          best_model <- model
          best_sig <- sig
          best_aic <- aic
        }
      }
    }
  }
  
  return(list(arch.model = arch.model, best_model = best_model))
}

model.arch1 <- arch.garch(data = train_nvda.diff, pp = 12, qq = 0, arma.order = arma.order)
datatable(model.arch1$arch.model, filter = 'top', 
          options = list(pageLength = 6))

cat("Best Model:", model.arch1$best_model)

## Best Model: ARCH(1)

model.arch2 <- arch.garch(data=train_nvda.diff, pp=12, qq=0, arma.order=arma.order,
                         mean=TRUE)
datatable(model.arch2$arch.model, filter = 'top', 
          options = list(pageLength = 6))

cat("Best Model:", model.arch2$best_model)

## Best Model: ARCH(1)

#Model dengan intersep
archSpec <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)),
  mean.model=list(armaOrder = arma.order, include.mean = TRUE) )
model.arch_nvda1 <- ugarchfit(spec=archSpec, data=train_nvda.diff, solver="hybrid")
coef(model.arch_nvda1)

##            mu           ma1         omega        alpha1 
## -5.087292e-05  4.308273e-02  7.454659e-07  2.989799e-01

#Model tanpa intersep
archSpec <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)),
  mean.model=list(armaOrder = arma.order, include.mean = FALSE) )
model.arch_nvda2 <- ugarchfit(spec=archSpec, data=train_nvda.diff, solver="hybrid")
coef(model.arch_nvda2)

##          ma1        omega       alpha1 
## 4.514680e-02 7.495716e-07 2.961780e-01

plot.ts(ts(model.arch_nvda1@fit$residuals), col=col.nvda[1], lty=1, 
        xlab="Periode (Hari)", ylab="Conditional Variance",
        lwd = 2, main="ARIMA(0,1,1)-ARCH(1) NVDA", cex.main = 1.5)

6.2 GARCH

Referensi : GARCH models with R programming : a practical example with TESLA stock

Ya walaupun diatas sudah dijelaskan “tidak usah melanjutkan ke garch”, tapi kami hanya ingin mencoba. “Apa yang terjadi jika kita tetap melakukan model GARCH?”. Model terbaik dari GARCH(p,q) ditentukan dari banyaknya paramter yang signifikan dan AIC nya yang terkecil.

model.garch <- arch.garch(data = train_nvda.diff, pp = 5, qq = 5, arma.order = arma.order)
datatable(model.garch$arch.model, filter = 'top', 
          options = list(pageLength = 6))

cat("Best Model:", model.garch$best_model)

## Best Model: GARCH(1,1)

Terlihat bahwa memang data tidak cocok untuk model GARCH(p,q) karena dari rentang lag 1 sampai 12 yang dihasilkan model terbaiknya adalah GARCH(12,12). Namun jika pada modelnya, hanya persentase signifikannya sangat kecil. Sehingga tidak mungin untuk dijadikan model terbaik.

6.3 Model Terbaik

Sehingga diperoleh model terbaik yakni ARIMA(0,1,1)+ARCH(1).

plot(model.arch_nvda1, which="all")

## 
## please wait...calculating quantiles...

Diperoleh mean model yaitu ARIMA(0,1,1) dan varian model yaitu ARCH(1) untuk digunakan dalam peralaman, pada uji dignostik model pada plot QQ-norm terlihat bahwa sebaran residual tidak mengikuti garis lurus sehingga dapat disimpulkan bahwa residual tidak menyebar normal.

Box.test((sisaan_nvda.da)^2, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  (sisaan_nvda.da)^2
## X-squared = 4.2891, df = 1, p-value = 0.03836

7.1 Plot Ramalan

Model ARIMA(0,1,1)

#Plot
plot(ramalan_nvda.da, main="NVDA", xlab="Periode(Hari)", cex.main=4, 
     lwd=2.5, col = col.nvda[1])

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa semua ramalan cenderung stabil hingga akhir periode.

Model ARIMA(0,1,1)+ARCH(1)

plot(ramalan_nvda.arch1, which= 1)

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa semua ramalan cenderung stabil hingga akhir periode.

#--Diff inv--
  #ARIMA
  ramal_nvda <- diffinv(data.ramalan_nvda.da, differences=1) +
    train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
  
  #ARCH + Intersep
  ramal_nvda.arch1 <- diffinv(data.ramal_nvda.arch1, differences=1) +
    train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
  
  #ARCH - Intersep
  ramal_nvda.arch2 <- diffinv(data.ramal_nvda.arch2, differences=1) +
    train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih

#--Transformasi Balik--
  #ARIMA
  ramal_nvda <- 1/(ramal_nvda)^2
  
  #ARCH + Intersep
  ramal_nvda.arch1 <- 1/(ramal_nvda.arch1)^2

  #ARCH - Intersep
  ramal_nvda.arch2 <- 1/(ramal_nvda.arch2)^2

#Data Frame
ramal <- data.frame(
  Date = data2 %>% dplyr::select(Date)  ,
  `ARIMA` = c(as.vector(train_nvda.ts), ramal_nvda[-1] ),
  `ARCH+Intersep` = c(as.vector(train_nvda.ts), ramal_nvda.arch1[-1] ),
  `ARCH-Intersep` = c(as.vector(train_nvda.ts), ramal_nvda.arch2[-1] ),
  Name = rep(c("NVDA"), 
              each = nrow(train_nvda) + nrow(test_nvda) )
)
colnames(ramal) <- c("Date", "ARIMA", "ARCH+Intersep", "ARCH-Intersep", "Name")

chart <- 
ggplot() +
# Akhir data
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data2$Date) )), 
             linetype = "dotted", color = "red") +
  #Label Data Asli 
  annotate( "rect", alpha=.1, fill="#4EC2C1",
            xmin=as.Date(min(data2$Date)), 
            xmax=as.Date(data2$Date[nrow(train_nvda)]),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Asli", size=10) + 
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(train_nvda)," Hari (", 
                   round(nrow(train_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) + 
  
  #Label Data Ramal
  annotate( "rect", alpha=.1, fill="violetred",
            xmin=as.Date(data2$Date[nrow(train_nvda)]), 
            xmax=as.Date(max(data2$Date)),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Ramal", size=10) +
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(test_nvda)," Hari (", 
                   round(nrow(test_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) +
  
  geom_vline( #Buat garis batas data
    xintercept = as.Date(data2$Date[nrow(train_nvda)]) , 
             linetype = "dotted", color = "black", linewidth = 2) +

#Time Series
  #Data RAMAL ARIMA
    #Dari data asli & ramal
  geom_line(data = ramal, linewidth=2,
            aes(x = Date, y = `ARIMA`, col = "ARIMA")) +
  #Data Ramal ARCH + Intersep
    #Dari data ramal saja
  geom_line(data = ramal[nrow(train_nvda):nrow(nvda),], linewidth=2,
            aes(x = Date, y = `ARCH+Intersep`, col = "ARCH+")) +
  #Data Ramal ARCH - Intersep
    #Dari data ramal saja
  geom_line(data = ramal[nrow(train_nvda):nrow(nvda),], linewidth=2,
            aes(x = Date, y = `ARCH-Intersep`, col = "ARCH-")) +
  #Data Test
  geom_line(data = test_nvda, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "Test")) +
  
  #Penanda
  geom_point(data = tail(train_nvda, 1), alpha = .5, 
             aes(x = Date, y = `Adj Close`), stroke=2,
             size = 15, shape = 21, color = "black", fill="violetred") +
  scale_colour_manual(values = c("ARIMA"=col.nvda[1],"ARCH+"="#539282",
                                "ARCH-"="#93C6B9","Test"=col.nvda[2]
                                 )) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Akurasi Ramalan Saham NVDA",
       subtitle = "Perbandingan Hasil Ramalan dengan data Asli\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(ramal$Date)), 
               as.Date(max(ramal$Date)) + 80)
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
  
#Bagi Data Train & Test
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)), 
                   y = max(nvda$`Adj Close`)* .58, 
                   yend = max(nvda$`Adj Close`)*.7), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)), 
                      y = max(nvda$`Adj Close`)*.7, 
                      label = max(train_nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +

#Data Uji
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date))+60, 
                   y = nvda$`Adj Close`[nrow(nvda)], 
                   yend = nvda$`Adj Close`[nrow(nvda)]), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date))+60, 
                      y = nvda$`Adj Close`[nrow(nvda)], 
                      label = "Data Asli" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +

# Data Ramal ARIMA
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)) + 60, 
                   y = max(ramal$`ARIMA`), 
                   yend = max(ramal$`ARIMA`)), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)) +60, 
                      y = max(ramal$`ARIMA`), 
                      label = "ARIMA" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  ) +
  
# Data Ramal ARCH + Intersep
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)) + 70, 
                   y = max(ramal$`ARCH+Intersep`), 
                   yend = max(ramal$`ARCH+Intersep`)), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)) +70, 
                      y = max(ramal$`ARCH+Intersep`), 
                      label = "ARCH dengan <br> Intersep" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +
  
# Data Ramal ARCH - Intersep
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)) + 70, 
                   y = ramal$`ARCH-Intersep`[nrow(nvda)], 
                   yend = ramal$`ARCH-Intersep`[nrow(nvda)]), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#DC4F26", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)) +70, 
                      y = ramal$`ARCH-Intersep`[nrow(nvda)], 
                      label = "ARCH tanpa <br> Intersep" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#DC4F26", box.color = "white", parse = TRUE
  ) +

#Tanggal akhir
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)), 
                   y = -Inf, 
                   yend = max(nvda$`Adj Close`)*.15), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)), 
                      y = max(nvda$`Adj Close`)*.15, 
                      label = max(data2$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  )

chart

#Export Chart
ggsave("06_Akurasi_model.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Dapat dilihat bahwa rata-rata harga saham Amazon diramalkan akan cenderung sedikit menurun setiap periodenya. Di sisi lain, saham NVIDIA diramalkan akan menaik pesat setiap periodenya. Sedangkan harga saham lainnya cenderung naik sedikit.

7.2 Perbandingan Data Aktual dengan Data Ramal

#Perbandingan
perbandingan_nvda.da <- matrix(data=c(head(test_nvda.ts, n=nrow(test_nvda)), 
                                      ramal_nvda[-1],
                                      ramal_nvda.arch1[-1],
                                      ramal_nvda.arch2[-1]),
                               nrow = nrow(test_nvda), ncol = 4)
colnames(perbandingan_nvda.da) <- c("Aktual","ARIMA", 
                                    "ARCH+Intersep","ARCH-Intersep")
datatable(perbandingan_nvda.da, filter = 'top', 
          options = list(pageLength = 5))

7.3 Akurasi

#Akurasi Model ARIMA
akurasi.nvda <- accuracy(ts(ramal_nvda[-1]), 
                         head(test_nvda.ts, n=nrow(test_nvda) )) %>%
  as.data.frame() %>%
  cbind(Model = ramalan_nvda.da[["method"]]) %>% 
  relocate(Model, .before = 1)
row.names(akurasi.nvda) <- "ARIMA"

#Akurasi Model ARCH + Intersep
akurasi.nvda.arch1 <- accuracy(ts(ramal_nvda.arch1[-1]), 
                         head(test_nvda.ts, n=nrow(test_nvda) )) %>%
  as.data.frame() %>%
  cbind(Model = model.arch1$best_model) %>% 
  relocate(Model, .before = 1)
row.names(akurasi.nvda.arch1) <- "ARCH+Intersep"

#Akurasi Model ARCH - Intersep
akurasi.nvda.arch2 <- accuracy(ts(ramal_nvda.arch2[-1]), 
                         head(test_nvda.ts, n=nrow(test_nvda) )) %>%
  as.data.frame() %>%
  cbind(Model = model.arch2$best_model) %>% 
  relocate(Model, .before = 1)
row.names(akurasi.nvda.arch2) <- "ARCH-Intersep"

akurasi <- rbind(akurasi.nvda, akurasi.nvda.arch1, akurasi.nvda.arch2) %>%
  mutate(Keterangan = case_when(
    MAPE < 10 ~ "Sangat Baik",
    MAPE >= 10 & MAPE <= 20 ~ "Baik",
    MAPE > 20 & MAPE <= 50 ~ "Layak",
    MAPE > 50 ~ "Tidak Akurat"
  )) %>% relocate(Keterangan, Model, MAPE)

datatable(akurasi, filter = 'top', 
          options = list(pageLength = 7))

Dari tabel diatas terlihat bahwa model ARIMA(0,1,1)-ARCH(1) dengan Intersep berkatergori “Baik” dengan nilai MAPE nya sebesar \(15.64\%\) (kurang dari \(20\%\)).

8.1 Stasioneritas

index <- seq(1:nrow(nvda)) #Beda titik potong
bc =boxcox(nvda.ts~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main = 2)

lambda_nvda.ramal <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.ramal <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK

lsk.ramal <- data.frame(
  Lambda = c(lambda_nvda.ramal ),
  Batas.Bawah = c( min(sk_nvda.ramal)  ),
  Batas.Atas = c(max(sk_nvda.ramal) )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk.ramal$Keterangan <- apply(lsk.ramal, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk.ramal) <- c("NVDA")

datatable(lsk.ramal, filter = "none")

tran.ramal <- lsk.ramal %>%
  mutate(
    Transformasi = trans(Lambda),
    `Transformasi Balik` = trans_balik(Lambda)
  ) %>% dplyr::select(Lambda, Transformasi, `Transformasi Balik`)

datatable(tran.ramal, filter = 'none')

#Transformasi
nvda.ts.new <- 1/sqrt(nvda.ts)

#Plot BoxCox
index <- seq(1:nrow(nvda))
bc <-boxcox(nvda.ts.new ~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main=2)

lambda_nvda.ramal.new <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.ramal.new <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK

lsk.ramal.new <- data.frame(
  Lambda = c(lambda_nvda.ramal.new ),
  Batas.Bawah = c( min(sk_nvda.ramal.new)  ),
  Batas.Atas = c( max(sk_nvda.ramal.new)  )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk.ramal.new$Keterangan <- apply(lsk.ramal.new, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk.ramal.new) <- c("NVDA")

datatable(lsk.ramal.new, filter = 'none')

#Diff
nvda.diff <- diff(nvda.ts.new, differences = 1) 

tseries::adf.test(nvda.diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  nvda.diff
## Dickey-Fuller = -8.9643, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

if (tseries::adf.test(nvda.diff)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
  cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}

## Karena p-value < 0.05, Maka Rataan Stasioner

8.2 Model

#Model
archSpec <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)), #ARCH(1)
  mean.model=list(armaOrder = arma.order, include.mean = TRUE) ) #ARIMA(0,1,1) + Intersep
model.ramal <- ugarchfit(spec=archSpec, data=nvda.diff, solver="hybrid")

#--Forecast--
nDate <- seq(from = max(data2$Date) + 1, 
             to = as.Date("2023-12-31"), #Dampe akhir tahun
             by = "days")

ramalan_nvda <- ugarchforecast(model.ramal, data = nvda.diff, 
                      n.ahead = length(nDate))
data.ramal_nvda <- ramalan_nvda@forecast[["seriesFor"]]

data.ramal_nvda.bb <- ramalan_nvda@forecast[["seriesFor"]]/
  -( ramalan_nvda@forecast[["sigmaFor"]] )^2

data.ramal_nvda.ba <- ramalan_nvda@forecast[["seriesFor"]]/
  ( ramalan_nvda@forecast[["sigmaFor"]] )^2
#--Diff inv--
  #Mean
  ramal_nvda <- diffinv(data.ramal_nvda, differences=1) +
    nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
  #BB
  ramal_nvda.bb <- diffinv(data.ramal_nvda.bb, differences=1) +
    nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
  #BA
  ramal_nvda.ba <- diffinv(data.ramal_nvda.ba, differences=1) +
    nvda.ts.new[length(nvda.ts)] #nilai akhir data latih


#--Transformasi Balik--
  ramal_nvda <- 1/(ramal_nvda)^2
  ramal_nvda.bb <- 1/(ramal_nvda.bb)^2
  ramal_nvda.ba <- 1/(ramal_nvda.ba)^2

#Data Frame
data_ramal <- data.frame(
  Date = seq(from = min(data2$Date),
             to = max(nDate), by="days")  ,
  `Adj Close` = c(as.vector(nvda.ts), ramal_nvda[-1] ),
  `BB` = c(as.vector(nvda.ts), ramal_nvda.bb[-1] ),
  `BA` = c(as.vector(nvda.ts), ramal_nvda.ba[-1] ),
  Name = rep(c("NVDA"), 
              each = length(nvda.ts) + length(nDate) )
)
colnames(data_ramal) <- c("Date", "Adj Close", "BB", "BA", "Name")

8.3 Plot

harga_akhir <- nvda$`Adj Close`[nrow(nvda)]
harga_ramal <- data_ramal$`Adj Close`[nrow(data_ramal)]
perc <- (harga_ramal-harga_akhir)/harga_akhir * 100

chart <- 
ggplot() +
# Akhir data
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data_ramal$Date) )), 
             linetype = "dotted", color = "red") +
  #Label Data Asli 
  annotate( "rect", alpha=.1, fill="#4EC2C1",
            xmin=as.Date("2022-06-01"), 
            xmax=as.Date(nvda$Date[nrow(nvda)]),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="#4EC2C1",
            x = as.Date(data_ramal$Date[ceiling(nrow(data_ramal)*.7)]), 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Saham Historis", size=10) + 
  annotate( "text", color="#4EC2C1",
            x = as.Date(data_ramal$Date[ceiling(nrow(data_ramal)*.7)]), 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(" "), size=10) + 
  
  #Label Data Ramal
  annotate( "rect", alpha=.1, fill="violetred",
            xmin=as.Date(data_ramal$Date[nrow(nvda)]), 
            xmax=as.Date(max(data_ramal$Date)),
            ymin=-Inf, ymax=Inf ) + 
  annotate( "text", color="violetred",
            x = as.Date(max(nDate)-25) , 
            y = min(data2$`Adj Close`) * 10, 
    label = "Data\nRamal", size=10) +
  
  geom_vline( #Buat garis batas data
    xintercept = as.Date(data_ramal$Date[nrow(nvda)]) , 
             linetype = "dotted", color = "black", linewidth = 2) +

#Time Series
  #Data RAMAL 
  geom_line(data = data_ramal, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "NVDA")) +
  
  #Penanda
  geom_point(data = tail(nvda, 1), alpha = .5, 
             aes(x = Date, y = `Adj Close`), stroke=2,
             size = 15, shape = 21, color = "black", fill="violetred") +
  scale_colour_manual(values = cols) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Ramalan Saham NVDA",
       subtitle = "Peramalan Saham NVDA Hingga Akhir Tahun 2023\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date("2022-06-01"), 
               as.Date(max(data_ramal$Date)) + 40)
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
  
#Tanggal akhir data asli
  geom_richtext(
    data = data.frame(x = as.Date(max(nvda$Date)), 
                      y = max(nvda$`Adj Close`)*.75, 
                      label = max(nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +
  
# Harga akhir Data Asli
  geom_segment(aes(x = as.Date(max(nvda$Date)), 
                   xend = as.Date(max(nvda$Date))- 30, 
                   y = nvda$`Adj Close`[nrow(nvda)], 
                   yend = max(data_ramal$`Adj Close`)*.9 ) , 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(nvda$Date)) - 30, 
                      y = max(data_ramal$`Adj Close`) *.9, 
                      label = paste0("$ ", harga_akhir %>% round(.,2)) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +
  
# Harga akhir Data Ramal
  geom_segment(aes(x = as.Date(max(data_ramal$Date)), 
                   xend = as.Date(max(data_ramal$Date)) + 30, 
                   y = max(data_ramal$`Adj Close`), 
                   yend = max(data_ramal$`Adj Close`)), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data_ramal$Date)) +30, 
                      y = max(data_ramal$`Adj Close`), 
                      label = paste0("$ ", harga_ramal %>% round(.,2)) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +
# Harga akhir data ramal beda brp persen
  annotate( "text", color=ifelse(perc>0, "green4", "red3"),
            x = as.Date(max(data_ramal$Date)) +30, 
            y = max(data_ramal$`Adj Close`) *.94, 
    label = paste0(ifelse(perc>0, "+ ", "- "), perc %>% round(.,2), "%"), 
    size=8 ) + 

#Tanggal akhir ramal
  geom_segment(aes(x = as.Date(max(data_ramal$Date)), 
                   xend = as.Date(max(data_ramal$Date)), 
                   y = -Inf, 
                   yend = max(nvda$`Adj Close`)*.15), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data_ramal$Date)), 
                      y = max(nvda$`Adj Close`)*.15, 
                      label = max(data_ramal$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  )

chart

#Export Chart
ggsave("07_Ramalan.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Model ARIMA(0,1,1)-ARCH(1) dengan intersep digunakan untuk meramalkan hingga periode akhir tahun 2023. Pada Gambar 12 terlihat bahwa adjusted close saham NVIDIA diramalkan akan menyentuh angka seharga 555,95 USD atau naik sebesar \(15,02\%\) hingga akhir tahun 2023.

9.1 Kesimpulan

Kesimpulan yang dapat ditarik pada penelitian ini adalah model ARIMA(0,1,1)-ARCH(1) dengan intersep merupakan model terbaik yang dapat meramalkan adjusted close saham NVIDIA dalam kurung waktu 11 November 2023 hingga 31 Desember 2023. Nilai MAPE yang didapatkan dari model terbaik adalah sebesar \(15.64\%\). Adjusted close saham NVIDIA diramalkan akan naik hingga menyentuh angka \(555.96\) USD atau naik sebesar \(15.02\%\). Hal ini menunjukkan bahwa jika ingin membeli saham harian NVIDIA pada kurung waktu 11 November 2023 hingga 31 Desember 2023 merupakan hal yang tepat dikarenakan saham akan naik.

9.2 Referensi

• Forecasting: Principles and Practice, Rob J Hyndman and George Athanasopoulos. Monash University, Australia.

• BOOSTEDML : Time series, Articles on Statistics and Machine Learning for Healthcare.